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Abstact: We study numerically and analytically the properties of the stationary state of a particle 
moving under the influence of an electric field E in a two dimensional periodic Lorentz gas with the 
energy kept constant by a Gaussian thermostat. Numerically the current appears to be a continuous 
function of E whose derivative varies very irregularly, possibly in a discontinuous manner. We 
argue for the non differ entibility of the current as a function of E utilizing a symbolic description 
of the dynamics based on the discontinuities of the collision map. The decay of correlations and 
the behavior of the diffusion constant are also investigated. 



Keywords: Thermostatted Lorentz gas, steady state current, smoothness, regularity, symbolic 
dynamics. 



1. Introduction 

There has been much interest in modeling stationary nonequilibrium states (SNS) of physical systems by closed 
dynamical systems evolving under a deterministic (time- reversible) non-Hamiltonian dynamics [H][EM][GC]. 
This is in contrast to modeling by "open" systems with some stochastic dynamics at the boundaries or by 
Hamiltonian coupling to infinite reservoirs [EPR1][EPR2][BL][LS]. Some such modeling is necessary to obtain 
SNS since the evolution of an isolated system will only permit equilibrium stationary states (or some very 
unstable measures). 

These dynamical system models modify the Hamiltonian time evolution by the addition of purely formal (i.e. 
having no connection with the actual dynamics of the system) thermostats such as a Gaussian thermostat. This 
keeps the kinetic or total energy of the nonequilibrium system, e.g. one subjected to external forces which do 
work on it and drive it away from equilibrium, constant in time thereby permitting a SNS to exist. Such a 
model of a stationary current carrying state produced by a steady electric field E acting on the charges in a 
conductor forming a circle (periodic boundary conditions were first introduced by Moran and Hoover [MH]). 
Various aspects of this model were later studied both analytically and numerically [LNRM] [CELS] . It is however 
still not clear how well this dynamical system with its nonphysical way of extracting the heat generated by the 
current really models the essential features of electrical conduction in a physical system. This seems relevant 
to deciding on the utility of this "thermostatted" approach in nonequilibrium statistical mechanics. To answer 
this question we present here and in subsequent works detailed numerical and analytical studies of a class of 
such models describing stationary states of current carrying systems [CELS]. The system studied in the present 
work is the Moran-Hoover [MH] model of a single particle moving among a periodic two dimensional array of 
fixed scatterers (Sinai billiard) in the presence of an external field E and a Gaussian thermostat which keeps its 
kinetic energy fixed, see Fig. 1. In [BDLR],[BDGL] we consider the generalization of this model to iV-particles, 
N > 2, see also [BGG]. 

The equations describing the motion of the particle, including elastic scattering with the obstacles, are: 
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Fig. 1: General billiard structure with scatterers of radius R± and R2 in a periodic box with side length I. 
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where F obs (q) represents (symbolically) the collisions with the obstacles, ( • ) is the particle equal to unity and 
E = Ee with |e| = 1. 

Eq.(l.l) and thus the full dynamics keeps the kinetic energy ^- fixed so that the motion takes place on the three 
dimensional surface of constant energy which is the direct product of a torus (with holes due to the obstacles) 
and the circle with |v| = vo- We can set v = vqv where v is a unit vector v — (cos??, sin The microcanonical 
or uniform distribution on this three dimensional energy surface is invariant under the dynamics when E = 0. 
For E ^ the dynamics is not Hamiltonian. There is in fact a phase space contraction 7 on the energy surface 
equal to (E • v) jv\ [MH] [GG] [CELS] . 

Defining a new dimensionless time ^P-t where / is a characteristic length of the system (say half the length of 
the torus) and scaling q by I, v by vo and E by -4 leave the dynamics unchanged. We shall therefore take from 

now on I — 2 and vq = 1 (we choose I = 2 because in this way we can see the obstacle as having centers at (0.0) 
and (1,1)). Then given any initial absolutely continuous density ipo($, q) we shall be interested in the time 
evolved density ^ t ($,q) in the limit when t — > 00 which describes the stationary noncquilibrium state (SNS) of 
this system. 

Moran and Hoover [MH] carried out extensive numerical studies of ipt , or more precisely of the induced measure 
v t on the Poincare cross section parametrized by the points P = (a, (3) e S 1 x [— it/2, n/2], a corresponding to 
the angle locating the point on the perimeter of the obstacle where the collision took place and (3 to the angle 
between the normal vector at a and the unit velocity vector v (see subsection 2.1 for a more precise definition). 
Following the discrete time trajectory P n (in their case there was only one obstacle per unit cell of the triangular 
lattice, in the case of the billiard in Fig. 1 one should add an integer component to P n , i.e. P n = (a, (3, s) where 
s e {0, 1} indicates which obstacle is hit) they found that the density of points for large n had a fractal 
appearance, i.e. the stationary measure appeared to be singular with respect to the uniform (microcanonical) 
measure (47r) _1 cos adad(j) on P which is stationary and approached (in a weak sense) as t — > 00 when E = 0. 

This problem was later studied analytically by Chernov et al. [CELS] who proved that starting with any initial 
density i>o(6,q) there does indeed exist a limiting measure v t — ► v\ for the induced measure on the Poincare 
section as t — > 00, whose Hausdorff dimension is, for E ^ 0, strictly less than two (the corresponding measure n\ 
on the energy surface has Hausdorff dimension less than three) . They also showed that the stationary current 
j(E) = (v) E = er E[l + o(|E|)] for E small enough |E| < E , E « 1 (possibly as small as 10~ 20 ) where (-) E 
is the average with respect to \i\ and ao (generally a tensor) is given by the Kubo formula. (It also follows 
from their analysis that j(E) is continuous for |E| < Eq.) Recent work by Wojtkowski [W] suggests that it may 
be possible to extend the results of [CELS] to larger value of |E|; he found a precise value E below wich the 
system is hyperbolic and above which it is not. 

In the present work we examine j(E) and other properties of the stationary measure v + for E G (0.025,2) 
and e = a;, i.e. the field is in the x direction and so j(E) = j(E)x. Carrying out a numerical evaluation of the 
Kubo formula for cr (i.e. an integral of the autocorrelation function at zero field) we find good agreement with 
the results of the simulation for j(E) as E tends towards zero. A similar investigation, with less precision for 
small values of E, was carried out for a triangular geometry in [LNRM], where the results, which look similar 
to ours for E not too small, were analyzed in terms of periodic orbits. Here we focus on the behavior of the 
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current for E £ (0.025, 0.5) where rather precise numerical simulations indicate nonsmooth behavior of j(E) vs. 
E. In particular we analyze the apparent singularities of as a function of E in terms of the change in the 
singularities of the map from P n to P n +i caused by grazing collisions as E changes. We also analyze the formal 
expression for the derivatives of j(E) obtained from the Kawasaki formula [CELS], and the equivalent Ruelle 
formalism [R]. We identify possible origins of singularities and argue for a function j(E) that is non differentiablc 
or, at most, admits a weak derivative with a dense set of discontinuities although the non rigorous nature of our 
argument does not permit a definitive statement (a similar statement should hold for the dependence of j(E) 
on e). 

It should be noted that, as shown in [ER] the stationary current is directly related to the sum of the stable 
and unstable Lyapunov exponents of the system via the equality (j(E) • E) = — [A S (E) + A"(E)]. There is also 
a relation between the Hausdorff dimension of n + and the Lyapunov exponents Dh(v + ) = 1 — A U (E)/A S (E) 
(remember that A s (0) + A"(0) = but A S (E) + A"(E) < for E ^ 0). 

In addition to j(E) we also investigate the diffusion constant d(E) relative to the drift as a function of E and 
compare it to dj(E)/dE, equality constituting in the limit E — > the Einstein relation proven in [CELS]. Finally 
we compute the decay of the velocity autocorrelation function for E ^ and analyze the numerical results using 
the methods developed in [GG] . We find that the decay continues to be exponential (this was proven for E = 
in [Y]) with an oscillatory behavior which changes qualitatively at or close to those values of E at which we see 
the nonsmooth behavior in the higher derivatives of j(E). 



2. Current versus field: results and discussion 

We first describe the results of our simulations and then discuss possible theoretical explanations of their most 
interesting features. 



2.1. Numerical results 

We used three different methods to compute the current as a function of the field for the model in Fig. 1 with 
i?i = 0.39, i?2 = 0.79 1 and |v| = 1. The first two are based on the fact that the system appears to be ergodic 
for E<2 (this has been proved only for the case when the field is very small, see [CELS]). This allows us to 
compute the current by choosing a "typical" initial point and evolving it for a very long time while measuring 
the time average of the instantaneous current. 

The simulations were done by computing a trajectory of 10 9 successive collisions. The only difficulty in the 
algorithm is in computing the successive collision times in a fast enough manner without missing collisions that 
are nearly tangent. Our algorithm is based on a Newton method: it takes around 20 hours to compute 10 9 
successive collisions on a i586 processor with a 350 Mhz clock under a Linux operating system. The values 
of k(E) = j(E)/E obtained by the simulation are represented by the filled points in Fig. 2. The error-bars 
are computed by running 10 independent initial conditions and looking at the maximum and minimum values 
obtained in this way. It is reasonable to believe that, while the value of the current averaged over a long time 
T goes to as the field goes to 0, the fluctuations about that average with respect to the initial point should 
be more or less independent of E. Considering the heavy numerical work needed to evolve 10 points for 10 9 
successive collisions we did it for E = 0.025, 0.05, 0.1, 0.125 and 0.175 and interpolated linearly for the error-bar 
at other points (the reason of the choice to leave out only E = 0.75, 0.150 and 0.2 is because we think nothing 
new happens there in a sense that will become clear in the following subsection). 

To obtain more confidence in our numerics we also computed k(E) using the Kawasaki formula which allows us 
to compute also re(0) (for a rigorous proof of this formula when E is very small see [CELS]). Letting Ss(t, X), 
with X = (q, v), be the evolution generated by (1.1) and calling J(X) the velocity component along the 
direction of E, then j(E) = (J(X))e with (-)e representing the average with respect to /ig, one has the 
Kawasaki formula[CELS], 



These values coincide with the one used in [BGG]. There a system with many particles was studied and the radia were set to R± — 0.4. 
Rl — 0.8 plus a particles' radius of r — 0.01. It is clear that if only one particle is present this is equivalent to a point particle moving 
among scatterers with radii Ri — 0.39 and R2 — 0.79. 
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Fig. 2: Conductivity as a function of the field E. The filled circles represent the results from the time average on a single long trajectory 
with their error-bar. The crosses represent the result from the Kawasaki simulation. No error-bar is reported in this case (except for the 
value at E — 0) to enhance readability. In any case the errors on the Kawasaki values are larger than the ones on the long time averages. 



where (-)o represents the average with respect to the microcanonical distribution corresponding to . Observe 
that eq.(2.1) reduces for E = to the well known Green-Kubo expression for the linear response to an external 
field. 

There are two problems in using eq.(2.1) numerically. The first is the necessity of a very good random number 
generator for the initial conditions and the second is the impossibility of integrating the correlation function 
for long times if one wants to have good statistics. We solved the first problem by using a R250 generator 
(see [NR]) and, to ensure a better independence, we decorrelated the initial points by evolving then for 40 
collisions without field. In this situation we observed that at a time equal to 50 times the mean free flight 
time the correlation function in eq.(2.1) is of the order of 10 -4 allowing us to say that the truncation error is 
probably much smaller than the statistical error due to the fact that we evolved only 2 • 10 6 initial conditions. 
Observe that in this case we can assume that our experiment consists in the independent sampling of a random 
variable so that we can estimate the statistical error by its standard deviation divided by the square root of the 
number of initial points. We observe finally that this method, although numerically demanding, is theoretically 
more reliable than the first method used to get the filled points in Fig. 2 due to the intrinsic instability of the 
dynamics under consideration, see [GG] for a more precise discussion. 

The results of this simulation are also plotted in Fig. 2 using crosses. We report the error bar only for the 
value at E — to maintain readability of the graphs. We note that the error-bars for E > are approximatively 
twice that at E = (simulations at E = arc much easier and faster). This means that all the values are in 
within the error bar of those computed with the previous method. 

Due to the symmetry of our system we clearly have that n{E) = n(—E) so, assuming that k(E) is a smooth 
function of E we should be able to fit it for small E by k(E) = k(0) + k"(0)E 2 . Such a fit is indeed possible 
(in particular one would find k(0) = 0.169 and k"(0) = 0.53) but it will pass only through the error-bars of the 
first three points. From the fourth point on the graph has a very linear appearance (to be precise we can fit it 
well by k(E) = 0.169 + 0.042.E). The problem with this linear fit is that it would produce a discontinuity of 
the first derivative of n{E) at some value of E £ (0.05,0.075). Moreover, due to the very large error-bars near 
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E = 0, it is easy to see that the linear fit for E > 0.075 passes through the error-bars of all the points leaving 
open the possibility of a discontinuity of the first derivative of j(E) at E = 0. 

Clearly the question of discontinuities in the derivatives cannot be decided on the basis of numerical simulations. 
Wc need an analytical argument (to be transformed hopefully into a proof) to decide this question. We present 
such an argument in the following subsection. But we first show in Fig. 3 a graph of the current for higher 
value of the field to see if other discontinuities of the derivative can be seen. We used the first method described 
before Fig. 2 to compute the current from E = 0.2 to E = 2.0 at steps of AE = 0.025. Each run consisted 
of 5 • 10 7 collisions and the error bars have been computed for values of E = 0.25 + 0.25i by running 10 initial 
condition and interpolated for the other points (when they are not visible it is because they are smaller than 
the points). 
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Fig. 3: Behavior of the current for higher value of the field. When no error-bar is visible it means that it is smaller than the symbol used. 



Based on Fig. 3 it seems consistent to assume that k(E) is continuous for E < 2. Its derivative on the contrary 
continues to change in a very irregular, possibly discontinuous, way. The first value of E at which a big change 
is seen is around E = 0.35. To better appreciate this we plot the conductivity for E = 0.2 to E = 0.5 in Fig. 
4. 

As can be easily seen the most probable site of discontinuous behavior of the derivative are at fields in the 
intervals (0.225,0.25), (0.325,0.35) and (0.375,0.425). After discussing the possible origin of non smooth be- 
havior of the expectation of a smooth function with respect to we will return to the behavior of j(E) at 
those values of E. 

We briefly report on a third method to compute the value of the conductivity that we used to be sure that the 
above results are not due to a bias in our Newton scheme. The idea is the same as the first method but instead 
of the fast Newton style algorithm we used a discrete time integrator (a fourth order Rungc-Kutta integrator). 
To avoid missing collisions that are almost tangent we had to choose a very small time step, i.e. St = 10~ 4 . 
The slowness of this algorithm permitted us to check only a few points and they were consistent with the other 
methods. 

We close this section observing that our simulations support the validity of a central limit theorem for our 
system. In fact if one looks at the behavior of the maximum and minimum values observed among the ten 
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Fig. 4: Current versus field for E £ [0.2,0.5]. The intervals (0.225, 0.25) and (0.325, 0.35) and (0.375, 0.425) possibly contain discountinuity 
of the first derivative. 

runs as a function of the running time T one finds that one can fit it with a function of the form C / y/T with 
C 2 ~ (J(Xf) E (J(X))%. 



2.2. Analytical discussion 

To study the properties of k(E) as a function of E we can use eq.(2.1). We know that the correlation functions 
in eq.(2.1) are uniformly integrable for E < E for some very small E , see [CELS]. This allows us to exchange 
limits, 



E^E' 



lim / dt(J(X)J(S E {t,X))} = 



/>OC 

^ dt(J(X) lim i J(S E (t,X))) = 



(2.2) 



f 

Jo 



dt(J(X)J(S E >(t,X))) 



and conclude that k(E) is a continuous function of E for E < Eq. Although the argument in [CELS] gives, as 
already stated, a very small value for E the numerical results, some of which are reported at the end of this 
subsection, support the validity of the above reasoning at least for E<1. 

Differentiating eq.(2.1) to obtain the higher derivatives of n{E) as a function of E is essentially equivalent to 
what is formally done in [R]. We obtain 

dn{E) '°° rt 



dt / dt'(J(X)F(S E (t',X))-Vs E{ t',x)J(S E (t-t',S E (t',X)))) . 



(2.3) 



dE jo jo 

where F(X) = x— (x ■ v)v is the derivative respect to E of the right hand side of eq.(l.l). This can be rewritten 
as: 
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= dt dt'(J(X)J(S E {t',X))J{S E {t,X)))) 
Jo Jo 



dn(E) 

dE (2.4) 

+ dt dt'(Vs E (t',x)J(X)-F(S E (t',X))J(S E (t,X)))) 
Jo Jo 

Observe that the second term on the right hand side can be rewritten as: 

rOO r-t 

dt dt'(\7 x J(X)A- 1 (t',X)F(S E (t',Xj)J(SE(t,X))))o = 
Jo Jo 

/•CO r-t 

\ dt \ dt'(xA- 1 (t',X)xJ(S E (t,X)))}o+ (2.5) 
Jo Jo 

I* CO r-t 

\ dt dt'(xA- 1 (t',X)-vJ(S E (t',X))J(S E (t,X))) 
Jo Jo 

where A(t,X) is the tangent flow generated by S(t,X), i.e. A(t,X) = dS g% X) ■ Eq. (2.4), together with 
analogous expressions for the higher order derivatives can be obtained easily from the formalism developed in 
[R]. Note that since the system has Lyapunov exponents \ U (E) > and one A s (E) < (at least for not too 
large E) one should expect that xA~ x (t' ', X)x ~ e - * 8 ^*. This makes the convergence of the integrals in (2.5) 
very problematic. We will return to these equations later after we consider this question from a different point 
of view 2 . 

It is convenient at this point to introduce a discrete time version of our dynamical system making precise what 
was already sketched in the introduction. For every point X in phase space we can define t(X) as the first 
time at which the particle collides with one of the scatterers. Let now S E (X) = S e (t(X) + , X) where the + 
indicates that we choose the velocity after the collision. It is clear that S E {X) restricted to the set T of points 
X = (q, v), such that q is on the surface of an obstacle and v is directed outward of the obstacle, defines a 
dynamical system (the set T is usually called a Poincare section). We can parametrize the collision points P 
in T by two angles and an integer, i.e. (a, (3, s) where s £ {0, 1} is if the point is on the obstacle with radius 
0.79 and 1 if it is on the obstacle of radius 0.39, a £ [— it,tt] is the angle between the direction of the field and 
the point on the scatterer at which the particle collides and (3 £ [— it/2, it/2] is the angle between the outgoing 
velocity and the normal direction to the scatterer, see Fig. 1. We can also write T = T° U T 1 where P £ T s iff 
P = (a, (3, s). With a slight abuse of notation we will still indicate our dynamical system in the new coordinates 
by S E (P) (note that for E = the invariant measure on T is proportional to cos (3df3da). 

As we already mentioned in the introduction the invariant measure /i^ for S E (t,X) will induce a measure 
for S{P) on T. The average of any function G(X) with respect to \i E can be obtained from v\ via: 



(G(X)) E = — / G(P)du+(P) (2.6) 

T E Jr 

where Q{X) = J^ X) G{S E {t,X))dt and r+ = J T r(P)du+{P). 

A very useful tool for the discussion of regularity properties of dynamical system in terms of external parameters 
consists in the construction of a conjugating map, i.e. we can look for a family of functions h EtE i{P) with the 
property that: 



S E {h E , E ,{P)) = h EiE ,(S E ,{P)). (2.7) 

see [BKL] for an example of such a construction for an Anosov systems. It is easy to see that for smooth 
dynamical systems (e.g. Anosov or axiom A) such a family exists, at least in a small neighborhood of any given 
E, and the average of a smooth function has the same regularity properties as the conjugation[JLl][BKL]. In 
our case the existence of a smooth family h E ^ E > of conjugation would not be enough to prove that if Q(P ) is a 
smooth function then (G)% = It^(^'^ u e(^') 1S a sm0 °th function of E. One would also have to show that the 
local unstable foliation as a function of E has good smoothness property 3 . We will not touch this second point 
and focus our attention on the conjugation, hoping that an answer to the problem of its existence will give us 



Much of the following material is the result of discussions that the authors had with G. Gallavotti, N. Chernov, L.-S. Young, D. Rucllc 
and C. Dcttman. We arc particularly indebted to G. Gallavotti for explaining to us a possible method to construct a Markov partition 
for the billiard, construction that forms the basis of our analysis of the conjugation map. 
Sec [BKL] for aprecisc discussion of this point. 
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enough information on the smoothness of the k(E). We will try to construct such a conjugation and show that 
this is impossible at least in the usual strong sense. The analysis of the problems that one encounter in this 
attemped construction will lead us to identify the possible origin of the observed behavior of k{E). 
The main problem with our system is that the map Se(X) is only piecewise smooth 4 . In fact Se(X) is smooth 
except at points whose image is tangent at collision 5 , i.e. the discontinuity set £ E of Se{P) is defined by the set 
of P such that Se(P) = (a 1 , ±7r/2, s) for some a'. It is clear that we can divide £ E into two parts £ E S = £ E DT S 
corresponding to the two obstacles. Moreover, due to the time reversibility of the dynamics it is easy to see 
that if (a, (3, s) g £\ then (a, — (3, s) is in the discontinuity set £ E X of S E X . The set ^q' 1 is shown in Fig. 5. 




This set is formed by a finite number of disjoint locally one dimensional manifolds with possibly a finite 
number of bifurcations (see Fig. 5). The bifurcation points are the points P where a multiple tangency 
happens, i.e. where the /3-component of both Se{P) and S E (P) is in {— ir/2, ir/2}. It is also easy to see that 
the boundary points of these manifolds are in/3 = 7r/2or/3 = — n/2 and that every branch of £ E can be 
represented by an increasing function a(/3). Although Fig. 5 has been generated with the help of the molecular 
dynamics program described in the previous section it would not be difficult to write analytical expressions for 
the manifolds in £ E . We note that from cq.(2.6) we get that (J(X)) E = {J{P))%/t e where J{P) is not a 
smooth function on T. This should not be a problem considering that J, as any function Q{P) obtained from a 
smooth function G{X) on phase space in the same way, is discontinuous exaclty on £ E and smooth everywhere 
else. 

It is also easy to observe that T E — £ E n £ E is a discrete set with a finite number I E of points Pi. We can 
assume that Ie is a piecewise constant function of E and one can write Pi = Pi (E) for every Pi £ T E with Pi (E) 
defined and smooth on every open interval on which I E is constant. Moreover the points Pi(E), together with 



Wc will use the term smooth in the vague sense of "sufficiently differentiablc" : at the level of rigor of the forthcoming discussion it is 
not worthwhile to specify exactly what kind of regularity wc need. This ambiguity should not undermine the undcrstandability of the 
following reasoning. 

The trajectories issuing from these points are usually called "grazing trajectories" 
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the bifurcation points, divide the manifolds to which they belong into connected segments and these segments 
form the sides of polygons (generically squares and triangles at this level) that form a partition of T. Again we 
can assume that, for E in a small enough interval, these segments and polygons can be smoothly parametrized. 
Fig. 6 presents a snapshot of this situation for E = 0. 



i 

E=Q 




-3.14 -1.57 0.00 1.57 3.14 a 

Fig. 6: Discontinuity set S^' 1 and Sq 1 ' 1 for the map Se(P) on the obstacle of radius 0.39 

It is clear that any conjugating map \ie.e 1 has to map £ E to £ E , and £ E } to £ E , . A possible first approximation 
of such an He,e' consists in defining h!" EE , (Pi(E))\ x = Pi(E'). It is clear that such a map can be extended 

to a smooth map h% ] EI : T ^ T such that h ( E ] E ,(P) = P + (E' - E)5h { E ] E ,(P) where 8h { E ] E , is some smooth 
function and h E ^ E ,(£ E ) = £ E , for a = ±1. 

Although the map h^ EE , is surely not a conjugation we can consider it as a starting point and reproduce the 
above construction using the sets of discontinuities of S E {P) and of S E 2 (P), £ E and £ E 2 respectively, obtaining 
a new map h E ^ E ,. Observe that £ E — £ E U £\ where P S £ E iff the /3-component of S E {P) is in {— 7r/2, 7t/2}. 
An example of the set £% U £ E 2 is given in Fig. 7 

It is well known that the sets £ E will become dense as N — > oo so one can hope that the above construction 
will create a sequence of functions h^ E E , eventually converging to a real conjugation when N — > oo. Moreover, 
considering that all the so constructed approximations are smooth one can hope that also the limit will be 
smooth. There is nevertheless a major problem that makes such a convergence impossible, at least in this first 
naive meaning, as we will see in what follows. 

Let us again consider our first approximation h^ EE ,. It is easy to realize that we can associate to each of the 

connected regions in which £ E S divides T s a pair of integer c = (n x , n y ) indicating the difference in coordinates 
between the center of the obstacle with which the particle collides and the one from which it starts. From 
now on we will see the dynamics as taking place on the universal covering of the torus and we will assume 
that the origin coincides with the center of one of the scatterers of radius 0.79. We will denote the set of this 
pair of integers by D E S and call it the set of possible "forward histories of lenght one" for reasons that will 
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I 

£=0 



1.570 




-3.14 -1.57 0.00 1.57 3.14 



Fig. 7: Discontinuity set S^' 1 and £ Q 2,1 for the map Se(P) on the obstacle of radius 0.395 



become clear in what follows. Observe that the map Se is continuous from one side at every point on £ E . This 
allows us to define the set D E S c D E S of the pair of integers associated with points in £ E s (these represent the 
obstacles that can be grazed starting from a given obstacle). In Fig. 8 we show, in the two upper pictures, the 
set D E for E = and E = 0.45. One can realize that two new obstacles can be grazed at a field E — 0.45 
whereas they were not at a field E = 0. The two lower pictures show the corresponding appearence of a new 
branch of £ E and the fact that a new connected component of T^£ E is created. The number on the figures are 
meant to help associate corresponding regions. The real modification of the set D 1 ^ 1 takes place at some value 
of E £ (0.375, 0.4) but we show the situation at E = 0.45 because it would otherwise be difficult to observe it 
due to the small size of the region denoted by 5 in the lower right picture in Fig. 8. It is interesting to note 
that this change takes place in one of the intervals where we located a possible discontinuity of the derivative of 
k(E). We believe that this phenomenon is at the origin of those discontinuities and we will argue in this sense 
in what follows. 

Let us call Ed the exact value of E at which the new possible histories appear. We can still construct our map 
h E l_ SE E d +&E-< f° r a ver y small SE, mapping the regions numbered from 1 to 4 in Fig. 5 for E = Ed — SE in to 
the corresponding ones for E = Ed + SE. This would create the first approximation h E ^_ SE E d +5E to wna t we 

can call a "partial conjugation" h Ed _$ E Ed+ $ E to be obtained as a limit of functions h^_ SE E d +5E constructed 
iterating the above scheme. Using this function to compute the average of a smooth function on T as a function 
of E would give a smooth result but we will commit an error due to the missed region corresponding to the 
new possible history (region 5 in Fig. 5). It is easy to realize that this region will have a size proportional to 
VSE and thus an area proportional to SE. We can assume, as a first approximation, that the error committed 
in neglecting such a region is proportional to this area which would give a discontinuity in the first derivative. 
We will now check if we can use this idea to account for the other points of discontinuity we have indicated in 
the previous section. First of all we observe that something very similar to what we discussed above happens 
in the interval (0.325,0.35) with the role of B X E played by D E °, as shown in Fig. 9. 

Although it seems to agree quite well with the numerical observations this picture is still very partial. In fact 
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Fig. 8: D 1 ^ 1 for E — and E — 0.45 with a detailed wiev of a part of the corresponding set E^ 1 showing the appearence of a new branch 
and a new connected component in T 1 \£^b1, 1 (numbers represent corresponding regions). Although the sizes are not to scale a larger point 
stands for an obstacle of radius R2. 




Fig. 9: D^ 1 for E = and E = 0.35. 



it can happen that the set X\ = E^ 1 changes structure although the sets S^ 1 remain structurally identical. 
It is possible to see from Fig. 6 that almost all the points in X\ are near bifurcation points for S^ 1 . Looking 
carefully at the picture (see Fig. 9) one realizes that there are indeed 4 intesection points that delimit a small 
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square. We observe that, as we did before, we can associate with every polygon in T bounded by £ E ' s U £ E ' S a 
sequence of length 2 (c_i, Ci) of pairs of integers Ci = {n xA , n y s) representing the center on which the particle 
collides forward and backward in time. We can call the set of all these pairs the possible "symmetric history of 
length two" and denote it by M E . If one follows the evolution of the above mentioned square as a function of 
the field one realizes that at a field in (0.05, 0.075) it disappears (see Fig. 9). We can again reason as before to 
deduce that this phenomenon will give rise to a discontinuity of the first derivative. The possible definition of 
the map h E ^ E , in this situation appears more complex and we will not discuss it here leaving this question and 
a formalization of the above picture to a forthcoming study. 




-0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 a -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 



Fig. 10: Detail of S^ 1 U £ E 1?1 showing the change of Ig although S^ 1 changes smoothly 

We observe here that the above considerations can be iterated. We can consider the set D E ' S formed by 
sequences of length N, (ci, . . . ,cjv) of pairs of integers representing the possible consecutive collisions that a 
trajectory experiences starting from the obstacle s (we can call it the set of possible "forward histories of length 
N"). With a reasoning analogous to the one above we can expect a discontinuity of the first derivative when 
one of this set changes. It is easy to see that, as N grows, the number of histories in the set D E increases 
exponentially and, due to the instability of the dynamics, new possible histories will appear or disappear very 
often. On the other hand each of the regions in which T is divided by the set £ E is exponentially small in N 
so we can expect them to create exponentially small discontinuities of the first derivative. A similar argument 
holds for the set I E and the set of symmetric histories of length 2N. This suggests that one should be able 
to define a "partial conjugation between E and E'" mapping T\Ee,e' to T\Ee>,e where T,e>,e and ^e,e> are 
two sets with area proportional to \E' — E\ and a very complex structure, possibly fractal. Considering our 
comments after eq.(2.7) this would suggests that n(E), as well as the average of any other smooth function of 
X, is not a smooth function of E, e.g. it has a dense set of discontinuities. If one accept the above picture on 
the structure of these discontinuities one can still conjecture that the function n(E) is Lipschitz continuous in 
E, i.e. \k(E) — k(E')\ < C\E — E'\, for a suitable constant C. This would imply that k(E) is differentiable 
almost everywhere with respect to the Lebesgue measure. 

We note here that in the interval (0.25, 0.275) several of the above phenomena take place (i.e. D E and D% 
change and a symmetric history of length two disappears) so that a precise discussion of what happens is not 
possible at this point. 

The above picture may seem unnatural and without hope of rigorous formalization and practical use. In this 
respect it is interesting to note that one can consider the polygons that are rectangles in Fig. 6 as part of a 
Markov (we use this notion in a much weaker sense than the usual meaning of a partition of T, in a possibly 
denumerable set of rectangles whose boundaries are mapped covariantly one to the other for a finite number 
of iterations). Indeed not all of them will have the necessary property of covariance but it is possible to give 
a theoretically simple algorithm to detect which can be considered as part of a Markov partition. The part of 
phase space not covered by these squares should be analyzed using the discontinuity set of S E . This will permit 
to obtain more squares and to cover a bigger portion of phase space. Clearly this algorithm can be iterated 
and, hopefully, will give rise to a complete Markov partition. In this situation our partial conjugation map is 
obtained by simply matching points with the same symbolic history. Although the above description seems 
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quite vague we think that the algorithm can be explicitly written and implemented on a computer. We hope 
to come back to this point in a forthcoming work. 

We conclude this section with two remarks: 

• One is typically used to think of conjugation as mapping the unstable manifold of one system into the unstable 
manifold of the other and similarly for the stable. Our proposed map, in some sense really does this. In fact 
the set £ E will get more and more parallel to the unstable local foliation of our system. An analogous comment 
holds for the Markov partition. 



d B (t) I 

| £=0.0, 0.025, J). 05, 0.075 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 



Fig. 11: Logarithm of orrelation function with respect to the the Lebesgue measure for fields E — (solid line), E — 0.025 (short dashed), 
E = 0.05 (long dashed) and E - 0.075 (dotted dashed). 



• An idea of this situation can come also from (2.5). Indeed it is reasonable to think that A E (t, X) will be given 
by a function that is smooth everywhere except on the discontinuity generated by the collisions that happen 
before time t. This implies that if X is part of a grazing trajectory for a given E$ we can expect A E (t, X) to be 
non smooth as a function of E at E = E a . Considering that the set of grazing trajectories varies smootly with 
E and that (2.5) contains an integral over the variable X we do not expect any non-smoothness of the current 
to take place before the set of grazing trajectories changes its structure. A dimensional analysis gives the same 
result as for the conjugation but we do not report it here because we do not believe that the integrals in (2.5) 
converge. We show instead in Fig. 11 the behavior of d E (t) — \og(\D E (t)\) with D E {t) = (J(X)J(S E (t, X))} E 
as a function of t and E. We see that all these functions have a maxima for t near 3 and near 4. The value 
at the maximum is a decreasing function of the field. In general the correlation function appears to be 
decreasing in E at fixed X except for the appearence of a new maximum around t = 3.5. We think that this 
new maximum can be connected with the appearence of a new symbol in D E and so with the discontinuity 
in the derivative of k(E) for small E. We leave this point as a speculation because checking it is beyond the 
scope of this paper. 
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3. Correlation function and diffusion coefficient 

In addition to D E {t) we studied the current-current correlation function in the steady state: 

C E (t) = (J(X)J(S E (t,X))) E - (J(X))% 

To compute this function we used a method similar to the one used for computing the Kawasaki formula. 
After having generated an initial point using the R250 generator we evolve it for 20 collisions with the field 
E = to better decorrelate the initial choice. We then switch on the field and let the system evolve for 50 
collisions in such a way that the final point can be considered as distributed according to the stationary SRB 
distribution \i E . We then use such a point to compute a trajectory segment of lenght 40 over wich we compute 
the correlation. For small values of the field the correlation function C E {t) appears very similar to the Kawasaki 
correlation function. Considering that we already plotted Co(t) = D (t) in Fig. 13 we plot the logarithm of the 
correlation function c E (t) = log(|C B (t)|) for E = 0.025, 0.05, 0.075 and 0.1. 




Fig. 12: Estimation of the decay of the correlation function for E — 0.025, 0.05, 0.075 and 0.1 using a linear fit for the maxima of c^(i) 



One of the most interesting question about such correlation functions is whether it decay exponentially. In the 
case E = it has been proven in [Y] and the proof can be extended for small (possibly very small) value of 
the field [C] . We use here a method developed in [GG] to analyze this question. In that paper the exponential 
decay rate was obtained by dividing the set of maxima of the function Co {t) into two groups and fitting them by 
a linear law. Due to the fact that our simulations are shorter than the one used in [GG] (mainly for the reason 
that computing the collision time for E ^ is much harder than for E = since the exact solution contains 
transcendental functions) we see two maxima only for one of the two groups used in [GG] : the maximum near 
t = 2 and the one near t — 4. While a fit for E = gives us a value consistent with the one obtained by [GG] 6 it 

Note that our results are reported in our adimensional time while in [GG] they were reported in unit of the mean free fly time 
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is evident that the slope of the fitting line becomes more negative (i.e. a faster exponential decay) when the 
field is switched on. A numerical value can be deduced from the fit reported in Fig. 12 but we do not plot them 
separately because the error bar on each point would be too large. 

As already noted before the structure of Ce(£) changes as E changes. Moreover it seems to us (although with 
our data we cannot really make a definitive statement) that one would still be able to divide the maxima of 
the correlation function into groups each of which well fitted by a linear law whose slope depends on E 7 . This 
behavior can be explained if one assumes that the correlation function is given by a sum of (quasi)-periodic 
functions modulated by decaying exponentials. Such a behavior is well established in the case of a diamond 
billiard by the numerical work in [ACG]. The hypothesis formulated at the end of the previous section will 
mean, in this setting, that the terms in the above sum are linked with the symbolic dynamics of the system. 
We hope to come back to this problem in a forthcoming work. 

We computed the correlation function Cs(i) also for values of the field in the interval (0.275,0.4). All the 
above discussion can be applied to this case in the same way and we do not show any graphs because no new 
information can be obtained from them. The integral of the correlation function Ce (t) gives us the xx-element 
of the diffusion matrix for the system, i.e. d{E) = J °° CE(t)dt. We conclude this paper reporting the value of 
this quantity obtained from the above data. 



K(E)\d(E) 
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k.(E) versus d(E) 



0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 



Fig. 13: d(E) versus k(E) showing good agreement near E — but clear non equality at E 7^ 



As for the correlation function we have the data for E £ (0.025,0.1) and E e (0.275,0.375). In Fig. 13 the 
value of the diffusion coefficient is plotted together with the value of k{E) obtained from the Kawasaki formula. 
It is clear that for small fields the two formulas give almost the same value. Such an equality for E ^ would 
support the possibility of the validity of a Green-Kubo relation out of equilibrium [GG][ES]. Our results show 
that the Kawasaki formula and the Green-Kubo relation give different results for E ^ 0. 



In [GG] all the two groups of maxima where fitted by linear laws with the same slope. It is unclear whether such a property can be 
expected when E ^ 
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4. Conclusions 



In this paper we presented results of numerical simulations as well as analytical evidence for a non smooth 
behavior of the current as a function of the field in a thermostatted single particle model of electrical conduction. 
We mainly studied the invariant measure v\ for the discrete time dynamical system obtained by considering 
as timing events the collisions of the particle with the obstacles. Our analysis is based on (at the present 
nonrigorous) construction of a conjugation map or rather a "partial conjugation map" between the dynamics 
at two different values of the field E. This construction strongly relies on the expectation that it is possible to 
analyze the dynamics using symbols directly connected to the discontinuities of the collision map. The picture 
can probably be substantiated by constructing in a explicit way a Markov partition for the dynamical system. 

Our analysis implies that the average with respect to of any smooth function on T or, more generally, the 
average of every smooth function on phase space with respect to [1^, is a not twice differentiable with respect to 
the electric field E. Moreover one can conjecture that the first derivative exists in a weak sense but is everywhere 
discontinuous. Such a behavior has already been observed in [KD] for a very simple model. 

While our analysis is not rigorous it has the advantage of being, at least in principle, contructive. The scheme 
for the construction of a "partial conjugation" or better of a Markov partition can lead to a rigorous proof of 
our assertions together with the possibility of conducting computer assisted experiments on the billiards. 

As noted earlier the recent formalism developed by Ruellc in [R] suggests the presence of the same kind of non 
smoothness but we do not see a possibility of using it in a convincing way to argue in one direction or the other. 
It seems to us that a form of conjugation would be in any case necessary to show the convergence of eq.(2.3). 

We also analyze the behavior of the velocity- velocity correlation function and of the diffusion coefficient. 
Applying the methods developed in [GG] we argue that the correlation functions decay exponentially also when 
the field is non zero. Moreover the rate of decay appears to increase with the field, at least for small values 
of E. The detailed structure of the correlation function undergoes interesting changes and we try to correlate 
these changes with the property of the invariant measure discussed above. Due to the numerical difficulties of 
computing correlation functions our discussion remains speculative. 

Finally we computed the diffusion coefficient. For E = the diffusion coefficient gives the linear conductivity. 
It has been argued that this relation could be valid also for a field different from 0, see [ES]. This would 
imply that the integral of the correlation function computed with respect to the Lebesgue measure, and 
with respect to the SRB measure fig, should be equal. Our results show that this is not the case so that no 
direct extention of a Green-Kubo formula to non-equilibrium is possible (see [G] for an interesting proposal of 
what exending Green-Kubo out of equilibrium means). Results similar to ours were already present in [ES]. 
Considering that no good reason has ever been given for the above equality to hold we consider this result as 
natural. 
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